Studies 1-4: Figures

Luhrmann, Weisman, et al.

In this notebook, we create figures for Studies 1-4.

source("./scripts_general/dependencies.R")
source("./scripts_general/custom_funs.R")
source("./scripts_general/var_recode_contrast.R")
source("./study1/scripts_s1/s1_var_groups.R")
source("./study2/scripts_s2/s2_var_groups.R")
source("./study3/scripts_s3/s3_var_groups.R")
source("./study4/scripts_s4/s4_var_groups.R")

setwd("./study1/analysis/")

The working directory was changed to /Users/anthrouser/Documents/Projects/Mind and Spirit overview paper/mindspiritquant/study1/analysis inside a notebook chunk. The working directory will be reset when the chunk is finished running. Use the knitr root.dir option in the setup chunk to change the working directory for notebook chunks.

source("../../scripts_general/data_load.R")

Parsed with column specification:
cols(
  .default = col_double(),
  date = col_date(format = ""),
  researcher = col_character(),
  country = col_character(),
  site = col_character(),
  religion = col_character(),
  subject_gender = col_character(),
  subject_job = col_character(),
  subject_schedule = col_character(),
  subject_livedhere = col_character(),
  subject_lang = col_character(),
  subject_marital = col_character(),
  subject_hs = col_character(),
  subject_liveswith = col_character(),
  servicesperweek = col_character(),
  study = col_character()
)
See spec(...) for full column specifications.
Parsed with column specification:
cols(
  .default = col_double(),
  date = col_date(format = ""),
  researcher = col_character(),
  country = col_character(),
  quad = col_character(),
  subject_gender = col_character(),
  subject_job = col_character(),
  subject_schedule = col_character(),
  subject_livedhere = col_character(),
  subject_lang = col_character(),
  subject_marital = col_character(),
  subject_hs = col_character(),
  subject_school = col_character(),
  subject_liveswith = col_character(),
  servicesperweek = col_character(),
  godexpviaawe_freq = col_character(),
  sleephabit = col_character(),
  researcher_date = col_date(format = ""),
  researcher_notes = col_character(),
  site = col_character(),
  religion = col_character()
)
See spec(...) for full column specifications.
Parsed with column specification:
cols(
  .default = col_double(),
  ctry = col_character(),
  wher = col_character(),
  recr = col_character(),
  whoc = col_character(),
  demo_chur = col_character(),
  demo_ethn = col_character(),
  demo_maj = col_character(),
  demo_pocc = col_character(),
  demo_rlgn = col_character(),
  demo_sex = col_character()
)
See spec(...) for full column specifications.
Missing column names filled in: 'X1' [1]Parsed with column specification:
cols(
  .default = col_character(),
  X1 = col_double(),
  epi_subj = col_double(),
  epi_demo_age = col_double(),
  epi_demo_ses_num = col_double(),
  epi_demo_howr_num = col_double(),
  epi_demo_affr_num = col_double(),
  epi_demo_tung_num = col_double(),
  epi_demo_ytng_num = col_double(),
  epi_demo_affr_cat = col_logical(),
  epi_demo_tung_cat = col_logical(),
  epi_demo_ytng_cat = col_logical(),
  epi_version = col_double(),
  epi_charc = col_double()
)
See spec(...) for full column specifications.
2 parsing failures.
 row       col expected  actual                           file
1025 epi_charc a double Unclear '../../study3/data/d_demo.csv'
1027 epi_charc a double Unclear '../../study3/data/d_demo.csv'
Missing column names filled in: 'X1' [1]Parsed with column specification:
cols(
  X1 = col_double(),
  epi_ctry = col_character(),
  epi_subj = col_double(),
  score = col_double()
)
Joining, by = c("epi_ctry", "epi_subj")
Missing column names filled in: 'X1' [1]Parsed with column specification:
cols(
  X1 = col_double(),
  epi_ctry = col_character(),
  epi_subj = col_double(),
  score = col_double()
)
Joining, by = c("epi_ctry", "epi_subj")
Missing column names filled in: 'X1' [1]Parsed with column specification:
cols(
  X1 = col_double(),
  epi_ctry = col_character(),
  epi_subj = col_double(),
  question = col_character(),
  response = col_double(),
  order = col_double(),
  question_text = col_character()
)
Missing column names filled in: 'X1' [1]Parsed with column specification:
cols(
  X1 = col_double(),
  epi_ctry = col_character(),
  epi_subj = col_double(),
  question = col_character(),
  response = col_double(),
  order = col_double(),
  question_text = col_character()
)
Joining, by = c("epi_ctry", "epi_subj", "question", "response", "order", "question_text")
Joining, by = c("epi_ctry", "epi_subj")
Column `epi_ctry` joining character vector and factor, coercing into character vectorJoining, by = "epi_subj"
Missing column names filled in: 'X1' [1]Parsed with column specification:
cols(
  .default = col_double(),
  p7_ctry = col_character(),
  p7_abs_check = col_character(),
  p7_dse_check = col_character(),
  p7_se_check = col_character(),
  p7_unev_check = col_character(),
  p7_exsen_check = col_character(),
  p7_por_check = col_character(),
  p7_mm_check = col_character(),
  p7_dem_sex = col_character(),
  p7_dem_pocc = col_character(),
  p7_dem_major = col_character(),
  p7_dem_ethnicity = col_character(),
  p7_dem_rur.urb = col_character(),
  p7_dem_affrd.basics = col_character(),
  p7_dem_religion = col_character(),
  p7_dem_church = col_character(),
  p7_dem_holy.tung.gif = col_character(),
  p7_abs_child.exp_cat = col_logical(),
  p7_abs_poetic_cat = col_logical(),
  p7_abs_tv.real_cat = col_logical()
  # ... with 162 more columns
)
See spec(...) for full column specifications.
Joining, by = c("p7_ctry", "p7_subj", "question", "response", "scale")
Joining, by = c("p7_ctry", "p7_subj", "question", "response", "scale")
Joining, by = c("p7_ctry", "p7_subj", "question", "response", "scale")
Joining, by = c("p7_ctry", "p7_subj", "question", "response", "scale")
Joining, by = c("p7_ctry", "p7_subj", "question", "response", "scale")
Joining, by = c("p7_ctry", "p7_subj", "question", "response", "scale")
Joining, by = c("p7_ctry", "p7_subj", "question", "response", "scale")
Joining, by = c("p7_ctry", "p7_subj", "question", "response", "scale")
Joining, by = c("study", "p7_ctry", "p7_subj", "abs_score", "cog_score", "ctl_score", "dse_score", "hall_score", "para_score", "por_score", "pv_score", "spev_score")
Joining, by = c("study", "p7_ctry", "p7_subj", "abs_score", "cog_score", "ctl_score", "dse_score", "hall_score", "para_score", "por_score", "pv_score", "spev_score")

Summary statistics

Spiritual Events scores

Study 1

d1 %>%
  ggplot(aes(x = country, y = spirit_score, color = country, fill = country)) +
  geom_jitter(height = 0, width = 0.25, alpha = 0.2, show.legend = F ) +
  geom_pointrange(data = . %>%
                    group_by(country) %>%
                    summarise(mean = mean(spirit_score, na.rm = T),
                              sd = sd(spirit_score, na.rm = T)) %>%
                    ungroup(),
                  aes(y = mean, ymin = mean - sd, ymax = mean + sd),
                  shape = 23, color = "black",
                  show.legend = F) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  scale_y_continuous(limits = c(0, 1), breaks = seq(0, 1, 0.25)) +
  labs(x = "Country", y = "Study 1: Spiritual Events score (range: 0-1)",
       caption = "Error bars are ±1 standard deviation from the mean")

Study 2

d2 %>%
  ggplot(aes(x = country, y = spev_score, color = country, fill = country)) +
  geom_jitter(height = 0, width = 0.25, alpha = 0.2, show.legend = F ) +
  geom_pointrange(data = . %>%
                    group_by(country) %>%
                    summarise(mean = mean(spev_score, na.rm = T),
                              sd = sd(spev_score, na.rm = T)) %>%
                    ungroup(),
                  aes(y = mean, ymin = mean - sd, ymax = mean + sd),
                  shape = 23, color = "black",
                  show.legend = F) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  scale_y_continuous(limits = c(NA, 4), breaks = seq(0, 4, 1)) +
  labs(x = "Country", y = "Study 2: Spiritual Events score (range: 0-4)",
       caption = "Error bars are ±1 standard deviation from the mean")

Study 3

d3 %>%
  # correct for scaling in original dataset
  mutate(spirit_score = spirit_score * 4) %>%
  ggplot(aes(x = epi_ctry, y = spirit_score, color = epi_ctry, fill = epi_ctry)) +
  geom_jitter(height = 0, width = 0.25, alpha = 0.2, show.legend = F ) +
  geom_pointrange(data = . %>%
                    group_by(epi_ctry) %>%
                    summarise(mean = mean(spirit_score, na.rm = T),
                              sd = sd(spirit_score, na.rm = T)) %>%
                    ungroup(),
                  aes(y = mean, ymin = mean - sd, ymax = mean + sd),
                  shape = 23, color = "black",
                  show.legend = F) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  scale_y_continuous(limits = c(NA, 4), breaks = seq(0, 4, 1)) +
  labs(x = "Country", y = "Study 3: Spiritual Events score (range: 0-4)",
       caption = "Error bars are ±1 standard deviation from the mean")

Study 4

d4 %>%
  ggplot(aes(x = p7_ctry, y = spev_score, color = p7_ctry, fill = p7_ctry)) +
  geom_jitter(height = 0, width = 0.25, alpha = 0.2, show.legend = F ) +
  geom_pointrange(data = . %>%
                    group_by(p7_ctry) %>%
                    summarise(mean = mean(spev_score, na.rm = T),
                              sd = sd(spev_score, na.rm = T)) %>%
                    ungroup(),
                  aes(y = mean, ymin = mean - sd, ymax = mean + sd),
                  shape = 23, color = "black",
                  show.legend = F) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  scale_y_continuous(limits = c(NA, 4), breaks = seq(0, 4, 1)) +
  labs(x = "Country", y = "Study 4: Spiritual Events score (range: 0-4)",
       caption = "Error bars are ±1 standard deviation from the mean")

All studies

d_all <- d1 %>%
  select(study, country, site, religion, subject_id, 
         por_score, abs_score, spirit_score) %>%
  rename(spev_score = spirit_score,
         pv_score = por_score) %>%
  mutate(religion = recode_factor(religion,
                       "indigenous" = "Indigenous Religion",
                       "charismatic" = "Charismatic Christianity"),
         # rescale to 0-1
         pv_score = pv_score/3) %>%
  full_join(d2 %>% 
              select(study, country, subj, 
                     abs_score, spev_score, dse_score) %>% 
              rename(subject_id = subj) %>%
              # rescale to 0-1
              mutate(spev_score = spev_score/4,
                     dse_score = dse_score/5,
                     religion = "General Population")) %>%
  full_join(d3 %>% 
              select(study, epi_ctry, epi_sample, epi_subj, 
                     por_score, spirit_score) %>%
              rename(country = epi_ctry, 
                     religion = epi_sample,
                     subject_id = epi_subj,
                     spev_score = spirit_score) %>%
              mutate(religion = recode_factor(religion,
                "general population" = "General Population",
                "charismatic" = "Charismatic Christianity"))) %>%
  full_join(d4 %>%
              select(study, p7_ctry, p7_subj, 
                     por_score, pv_score, abs_score, spev_score, dse_score) %>%
              rename(country = p7_ctry,
                     subject_id = p7_subj) %>%
              # rescale to 0-1
              mutate(por_score = por_score/2,
                     pv_score = pv_score/3,
                     spev_score = spev_score/4,
                     dse_score = dse_score/5,
                     religion = "General Population")) %>%
  mutate(religion = factor(religion, 
                           levels = c("General Population", 
                                      "Indigenous Religion",
                                      "Charismatic Christianity")),
         study = gsub("study", "Study", study)) 

Joining, by = c("study", "country", "religion", "subject_id", "abs_score", "spev_score")
Column `religion` joining factor and character vector, coercing into character vectorJoining, by = c("study", "country", "religion", "subject_id", "spev_score")
Column `religion` joining character vector and factor, coercing into character vectorJoining, by = c("study", "country", "religion", "subject_id", "pv_score", "abs_score", "spev_score", "dse_score", "por_score")

d_all %>%
  ggplot(aes(x = country, y = spev_score, color = country, fill = country,
             group = study)) +
  geom_point(position = position_jitterdodge(jitter.width = 0.25,
                                   jitter.height = 0,
                                   dodge.width = 0.75), 
             alpha = 0.1) +
  geom_pointrange(data = d_spev %>%
                    group_by(study, country) %>%
                    summarise(mean = mean(spev_score, na.rm = T),
                              sd = sd(spev_score, na.rm = T)) %>%
                    ungroup(),
                  aes(y = mean, ymin = mean - sd, ymax = mean + sd, shape = study), 
                  position = position_dodge(width = 0.75),
                  color = "black") +
  scale_shape_manual(values = 21:24) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  guides(color = F, fill = F, 
         shape = guide_legend(override.aes = list(fill = "black"))) +
  theme(legend.position = "bottom") +
  labs(x = "Country", y = "Spiritual Events score (rescaled to 0-1)",
       shape = "Study", 
       caption = "Error bars are ±1 standard deviation from the mean")

d_spev_sum <- d_all %>%
  filter(!is.na(spev_score)) %>%
  group_by(study, country, religion) %>%
  summarise(mean = mean(spev_score, na.rm = T),
            sd = sd(spev_score, na.rm = T),
            n = n()) %>%
  ungroup()

d_all %>%
  ggplot(aes(x = country, y = spev_score, 
             color = country, fill = country,
             group = religion)) +
  facet_grid(~ study, scales = "free", space = "free") +
  geom_point(position = position_jitterdodge(jitter.width = 0.8,
                                    jitter.height = 0.02,
                                    dodge.width = 0.75), 
             alpha = 0.15) +
  geom_pointrange(data = d_spev_sum,
                  aes(y = mean, ymin = mean - sd, ymax = mean + sd, 
                      shape = religion), 
                  position = position_dodge(width = 0.75),
                  fill = "black",
                  color = "black") +
  geom_text(data = d_spev_sum %>%
              mutate(ypos = case_when(
                grepl("charismatic", tolower(religion)) ~ mean + sd + 0.05,
                TRUE ~ mean - sd - 0.05)),
            aes(y = ypos, label = paste0("n=", n)), 
            position = position_dodge(width = 0.75),
            size = 3, color = "black") +
  scale_shape_manual(values = 21:24) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  guides(color = F, fill = F, 
         shape = guide_legend(override.aes = list(fill = "black"))) +
  theme(legend.position = "bottom",
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
  labs(x = "Country", y = "Spiritual Events score (rescaled to 0-1)",
       # caption = "Error bars are ±1 standard deviation from the mean",
       shape = "Religion")

Relationships

All studies

d_all %>%
  gather(spirit_scale, spirit_score, c(spev_score, dse_score)) %>%
  gather(poros_scale, poros_score, c(por_score, pv_score)) %>%
  mutate(spirit_scale = recode_factor(spirit_scale,
                                      "spev_score" = "Spiritual Events",
                                      "dse_score" = "Daily Spiritual Experiences"),
         poros_scale = recode_factor(poros_scale,
                                     "pv_score" = "Porosity Vignettes",
                                     "por_score" = "Porosity Scale"),
         study_scale = paste(study, poros_scale, sep = ": ")) %>%
  filter(!is.na(poros_score)) %>%
  ggplot(aes(x = poros_score, y = spirit_score)) +
  facet_grid(spirit_scale ~ study_scale) +
  geom_point(data = . %>% distinct(study, study_scale, country, 
                            poros_scale, poros_score,
                            spirit_scale, spirit_score),
             aes(color = country), alpha = 0.1) +
  geom_smooth(aes(color = country), method = "lm", 
              lty = 2, size = 0.7, alpha = 0, show.legend = F) +
  geom_smooth(method = "lm", color = "black", alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  theme(legend.position = "bottom", 
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
  guides(color = guide_legend(override.aes = list(alpha = 1))) +
  labs(x = "Score on porosity measure (rescaled to 0-1)",
       y = "Score on spiritual experience measure (rescaled to 0-1)",
       # caption = "Solid black lines correspond to to the overall trend, collapsing across countries",
       color = "Country")

d_all %>%
  gather(spirit_scale, spirit_score, c(spev_score, dse_score)) %>%
  mutate(spirit_scale = recode_factor(spirit_scale,
                                      "spev_score" = "Spiritual Events",
                                      "dse_score" = "Daily Spiritual Experiences"),
         study_scale = paste(study, "Absorption scale", sep = ": ")) %>%
  filter(!is.na(abs_score)) %>%
  ggplot(aes(x = abs_score, y = spirit_score)) +
  facet_grid(spirit_scale ~ study_scale) +
  geom_point(data = . %>% distinct(study, study_scale, country, 
                            abs_score,
                            spirit_scale, spirit_score),
             aes(color = country), alpha = 0.2) +
  geom_smooth(aes(color = country), method = "lm", 
              lty = 2, size = 0.7, alpha = 0, show.legend = F) +
  geom_smooth(method = "lm", color = "black", alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  theme(legend.position = "bottom", 
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
  guides(color = guide_legend(override.aes = list(alpha = 1))) +
  labs(x = "Score on absorption measure (rescaled to 0-1)",
       y = "Score on spiritual experience measure (rescaled to 0-1)",
       # caption = "Solid black lines correspond to to the overall trend, collapsing across countries",
       color = "Country")